Fast construction of the Kohn-Sham response function for molecules 
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Abstract 

The use of the LCAO (Linear Combination of Atomic Orbitals) method for excited states involves products of 
orbitals that are known to be linearly dependent. We identify a basis in the space of orbital products that is local for 
orbitals of finite support and with a residual error that vanishes exponentially with its dimension. As an application of 
our previously reported technique we compute the Kohn-Sham density response function xo for a molecule consisting 
of N atoms in N 2 N LJ operations, with N w the number of frequency points. We test our construction of xo by computing 
molecular spectra directly from the equations of Petersilka-Gossmann-Gross in N 2 N^ operations rather than from 
Casida's equations which takes N :i operations. We consider the good agreement with previously calculated molecular 
spectra as a validation of our construction of xo- Ongoing work indicates that our method is well suited for the 
computation of the GW self-energy E = iGW and we expect it to be useful in the analysis of exitonic effects in 
molecules. 

Accepted for publication in Physica Status Solidi, 05.02.2009 both as contribution to the Proceedings of TNT2009 
conference (Barcelona, Spain) and as a featured article. 

1 Introduction 

The method of "linear combination of atomic orbitals" (LCAO) goes back to the early days of quantum mechanics [1] 
and remains a good choice in ab-initio calculations of molecules or solids. LCAO provides a parsimonious basis for 
representing molecular orbitals and one-particle Green's functions. However, electronic structure theory also involves 
the electronic density and the expression for the electronic density contains all non vanishing products of orbitals, a 
set of quantities that are known to be linearly dependent. 

The many existing constructions of a basis in the space of products may be divided into two classes. A first 
type of construction, called "resolution of identity" method was proposed by Boys and Shavitt j2] and focuses on the 
representation of the electronic interaction, see also Casida [3|. For Slater-type functions a "density fitting" procedure 
was developed by Baerends el al |H [5] . The goal of both of these methods is to represent the electronic density (a 
two-center quantity within LCAO method) in a basis of one-center functions. For a good discussion of these techniques, 
see [6]. 

A second class of methods provides a basis for both response functions and interactions. For solids and within 
the muffin-tin approach, a "product basis" was proposed by Aryasetiawan and Gunnarsson [7] who removed the linear 
dependence by orthogonalizing an overlap matrix. Gaussian-type auxiliary functions were also used in solid state theory 
[5]. Blase and Ordejon [§] used Gram-Schmidt orthogonalization to eliminate the linear dependence from products on 
the same atom. 

Our own construction [101 111] belongs to the second class of methods. It was developed in the context of LCAO 
for orbitals of finite support and keeps the locality properties of the underlying atomic orbitals. In the present paper 
we further test our calculational framework by computing spectra of medium sized molecules and we also give an 
expression of the GW self-energy in our product basis. 

This paper is organized as follows. In the next section, we describe our construction of a basis in the space of 
products in terms of "dominant products" . A construction of the Kohn-Sham response function is described in section 
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Figure 1: On the left panel, the orbital products (of total azimuthal angular momentum m z — 0) are shown along the 
line connecting two (carbon) atoms. On the right panel, the corresponding dominant products are shown. The origin of 
x-axis is the middle point between atoms (on both panels). The linear dependence of the original orbital products and 
the small size of the dominant products basis are clearly visible. 



[3] and the application of this response function to the computation of molecular spectra is given in section [4] Finally, 
in section [5] we derive an expressions for the self-energy in Hedin's GW approximation in our basis and give our 
conclusions. 



2 Reducing the number of orbital products 

Restated in the succinct language of second quantization, the LCAO method consists of an expansion of the electron 
creation and annihilation operators if> + (r,t), ip(r,t) in terms of electron operators c a {t) that belong to atomic orbitals 

V + (r,t)~]Tr(r) C +(t). (1) 

a 

The functions {f a (r)} are atom centered orbitals, the operator c„ (i) create localized electrons in these orbitals and the 
symbol ~ alludes to the question of completeness of these LCAO orbitals, a difficulty we shall ignore in the following. 
Most molecular response functions involve the electronic density n(r,t) and, therefore, the square of the previous 
expansion 

n(r,t) = i>+(r,t)^(r,t) ~ £ f a (r)f b (r)d(t)c b (t). 

a,b 

It is well known in quantum chemistry that the products of orbitals {f a (f)f b (r)} that parametrize this density are 
linearly dependent. A good illustration of this fact is provided [T^] by the harmonic oscillator and its Hermite wave 
functions (f> n (x) = H n (x)e~ x ^ 2 where iv< - J ^ +1 ^ products {$ m {x)<t>n{x)}m,n=i..N span a space of only 2iV + l dimensions. 
The set of products {f a (r)f b (r)} is usually parametrized by extra auxiliary fitting functions j4j[5[|T3]. In this paper, 
we use an alternative procedure that involves no fitting functions whatsoever [TO]. We proceed in two steps: 

• for each pair of atoms, the orbitals of which overlap, we enumerate all products F M (r) = f a (r)f b (r), 

• we compute their metric or matrix of overlaps G MN , diagonalize this matrix and employ its eigenvectors and 
eigenvalues to define dominant products F x (r) 

G MN = f^P^P dW , (2) 
J \r — r'\ 

G MN X^ = \X M ; F x {v) = X x M F M {r). 

We use the Coulomb metric, the favorable properties of which are well known in quantum chemistry |14] , Although 
the procedure is carried out separately for each pair of atoms at a time, the overall result is a basis in the space of 
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all products. Our procedure provides us with (i) a set of dominant functions {F x (r)} and (ii) their relation with the 
original products 

/°(r)/ 6 (r)~ £ V x ab F x (r). (3) 

A>A min 

Here we ignore the products that have a Coulomb norm less than A m i n . The finite support of the original atomic 
orbitals and the locality of our construction are reflected in the sparse character of the "vertex" V£ . 

Due to (i) the finite support of the original orbitals and (ii) the locality of our procedure, identifying the dominant 
functions in a molecule of N atoms takes asymptotically only O(N) operations. Moreover, a very accurate represen- 
tation of orbital products as an expansion about intermediate points of each pair is possible thanks to two powerful 
algorithms developed by Talman [15] . For the case of bilocal products, the results of our procedure are illustrated in 
figure [l] 

For reasons not yet entirely understood [TH], the eigenvalues of the metric Q are asymptotically evenly spaced on 
a logarithmic scale, in formal analogy with the 1// noise that occurs in electronics [17] • As a welcome consequence, 
the residual (asymptotic) error of our algorithm vanishes exponentially fast with the number of dominant products 
retained. 

3 Computation of the Kohn— Sham density response in O(N 2 N L0 ) 
operations 

The usual LCAO method provides us with a tensor basis for the electron operators c a (t), Cj (t') and for their Green's 
function iG a f>(i — t') — T{c a (i)c£ (t')) . Similarly, our dominant products provide a tensor basis for the density 

A 

with n x {t) = ^c+(t)V?*<%(t). 

a.b 

Recalling that the response function coincides with the density-density correlator |18l 119] 

Xo(r,r',t-t') = -i{T{n(r,t)n(r' ,t')}) conncctcd , 
we obtain a representation of the response function 

in terms of a time-dependent matrix X«^(^ — 

X%(t - = -i(T{n M (t)n„(t')}>coBnected = i Y, v ^Gac{t - t')V: d G db {t' - t). (4) 

a,b,c,d 

From figure pi and from the locality of the vertex we see that for iV 3> 1 atoms, the computation of y°^(t — t') 




Figure 2: Particle-hole diagram for Kohn-Sham response in a basis of dominant products. The vertex V® b connects pairs 
of orbitals a, b to a dominant product /x. The propagators connect the orbitals within the atom quadruplet. For a given 
pair of dominant products (p,v), only the orbitals that belong to the corresponding atom quadruplet must be summed 
over. Therefore, the total computational effort scales as OlN 2 ^). 
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takes 0(N 2 ) operations. More precisely, the number of operations is 0{N 2 N u} ]ag(N uj y\ for N u frequencies where the 
Nu log(7V aj ) factor is due to the use of fast Fourier techniques in evaluating equation Q. 

To achieve a well controlled calculation, we express x%v ( w ) m terms of its spectral representation 

where a M „(A) is a spectral function associated with response function X°i/( w )' We now indicate how expression ^ can 
be derived, for instance, by representing the Green's functions in equation Q in terms of its spectral functions. We 
first recall the standard expression [181 120] for Green's function in terms of molecular orbitals ipE(r) 



G(r, r',t - t') = i${t' -t)J2 ipE{r)ipE{r')e- iE{t - t ' ) - iS(t - t') ^ i) E (r)^ E {r')<r 11 



LE(t-t') 

v^y jyE\> 

B<0 B>0 



We then rewrite it using the LCAO expression for the molecular orbitals tpE(r) = J3 X s f a (r) and introducing density 
matrices of particles p+ (s) = J2 E >o X a^5{s - E) and holes p~ b (s) = J2 E <o X aX^S(s - E) 

G ab (t - = m - t') J ds pMe-^'-V - i9(t - t') J ds p+We-W*-*). 
Inserting the last expression into equation Q, we identify the spectral function a M „(A) in equation (JsJ 

o^(a)= {vt a pi c ®v; d p M }(\). (6) 

a,b,c,d 

Here the overlined density ~p refers to a reflection of the argument p(x) = p{—x) and the notation <g> represents the 
convolution. 

The convolutions in equation |6| can be done in 0(N UJ log(N UJ )) operations using the FFT technique, provided the 
density matrices /oi(s) are discretised on an uniform grid. The discretisation of stick-like density matrices is done with 
a simple algorithm — dividing the spectral weight between adjacent grid points according to the distance between an 
eigenenergy E and the grid point [11] . The integration over A in equation (|5| can also be represented as a convolution 
and may be computed in a fast manner using the FFT technique [llj . 

We checked our results against the exact but slow expression of the Kohn-Sham density response, i. e. 



u) — (E — F) — \e(riE — uf) 



and its corresponding expression in our tensor basis (see figure [3] for a comparison between exact and fast results). This 
summarizes our method of constructing the Kohn-Sham density response function xo in 0(N 2 N UJ log^^)) operations. 
Next we will test our construction on molecular spectra. 

4 Application of xo to molecular spectra 

The most direct application of \o as constructed in the previous section is the computation of excitation spectra of 
molecules. 

The excitation properties of molecules are determined by an interacting response function \ that is related to the 
non interacting response function xo [22] via a Dyson type equation 

X = ~FT7 = 1 f — Xo- (7) 

OVcxt 1 — /HxcXO 

The TDDFT kernel /h xc = — ? becomes a matrix in the basis of dominant products 

on 

f " KC = JJ \r-r'\ J J W (r) (8) 

In this work, we focus on the local density approximation (LDA) for the exchange-correlation potential. In this case 
/xc(f, **') = dV d^"^ {f)S(r — r') and its matrix elements are local in coordinate space. 
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Figure 3: An element of the response function x^y of benzene. The Kohn-Sham eigenstates were generated using the 
SIESTA package [21] with default settings. The discretised response function is computed in the frequency window 
lj < w max , w max = 2 Rydberg with N u — 512 data points, e is chosen to be 1.5Ao;, where the discretisation spacing is 
Au> = 2w roax /iV w . 



4.1 Calculation of the interaction kernel fn xc 

In the construction of dominant products, we distinguish between (i) coincident and (ii) distinct or bilocal pairs of 
atoms. The former have full rotational symmetry, while the latter have only axial symmetry with respect to a line 
connecting their centers. Therefore the expansion of the bilocal products is done in an appropriately rotated coordinate 
frame. Both local and bilocal dominant products are expressed as a sum over angular-radial functions in a suitable 
coordinate system 

F"(r) =Y, F t(\r - C,\)Y lm ^(R,(r - C M )), (9) 

i 

where the rotation Ji M and the shift C M are determined by the atom pair. Using the theory of angular momentum, one 
can reduce the Coulomb interaction JJ — fejj^rp- - d 3 rd 3 r' to a sum over elementary two-center Coulomb integrals 

(C\C) = Sim(T " C r )& ';; | (T ' , ' C,) d 3 rdV (10) 

between two functions of spherical symmetry g; m (r). When the orbitals overlap, the Coulomb integrals are computed 
in the momentum space (where they become local) while, if they do not overlap, the Coulomb interaction is calculated 
exactly in terms of moments. For converting between coordinate and momentum space, we use Talman's fast Bessel 
transform |15j . 

Unlike the Hartree kernel, the LDA exchange-correlation kernel is local in coordinate space and a numerical 
integration is an appropriate procedure for it. The integrand F Al (r)/ xc (r)i ? "(r) is of non trivial support due to the 
lens-shaped support of the dominant products. However, the non spherical part of this volume is small for neighboring 
atoms and we used spherical coordinates centered about a midpoint between the centers of the dominant products 
F M (r), F v (r). The integration is done with a Gauss-Legendre method along the radial coordinate and with a Lebedev 
quadrature [23] over the solid angle. A moderate number of integration points leads to a sufficient accuracy in the 
exchange-correlation matrix elements. 

The computational complexity of the Hartree kernel and the exchange-correlation kernel fj£? are 0(N 2 ) and 
0(N), respectively. In practice, the calculation of the kernel /h xc is faster than the calculation of the non interacting 
response xo by an order of magnitude. 

4.2 An iterative method for finding the interacting polarizability 

Expression (TtJ) for the interacting response Xv-v involves matrix multiplications and inversion and this would appear to 
require 0(N^NJ) operations, more than the 0(N 2 N LJ ) complexity of the non interacting response X%,w Fortunately, 
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the polarizability P is an average quantitjj^] 

P{u>) =Y,d»X^)d u , (11) 

that is easy to compute using an iterative method of the Krylov type |TT] . A biorthogonal Lanczos method 24 allows 
us to find a simple representation for the (non hermitian) matrix A = 1 — /hxcX° 

A = Y,\L n )tnm{R m \, (12) 

where the matrix t nm is tridiagonal. The set of vectors \L n ) and (7? m | build up the identity in the Krylov space, 
{R m \L n } — 5 mn . Therefore, if we choose the starting vectors IL 1 ) and {R 1 ] in the direction of the dipole moment d 
and of xod, respectively 

\L l ) =d; (R'\ = X °d, (13) 
then the interacting polarizability takes a particularly simple form 

PH = ^1x0 MM). (14) 

Here, we multiply the Kohn-Sham (non interacting) polarizability with the first element of the inverse matrix t nm . 

The computational complexity of the iterative method is 0(N 2 Nk), where Nk is the dimension of the Krylov 
subspace. In our calculation Nk is very small compared to the number of dominant products and grows slowly as we 
increase the size of the system. In conclusion, our method requires a total of 0{N 2 N^ log(A^)) operations to compute 
the molecular polarizability P(lo). 



4.3 Memory requirements 

The Kohn-Sham response matrix requires 0{N 2 N UJ ) memory. For instance, in the case of indigo dye, we need approxi- 
mately 4 GBytes, where we have exploited the diagonal symmetry of the response matrix, using single precision complex 
numbers, and calculating for 256 frequencies. This is beyond the capacity of most contemporary desktop computers. 
The storage of the response matrix on a hard disk is not an option because the matrix is generated element by element 
for all frequencies at once, but needed (in the iterative computation of polarizability) as a frequency dependent matrix. 
For this reason, we have to employ parallel machines if a large molecule, say 100 atoms, is to be treated. Currently, 
we are implementing OpenMP and MPI parallelized versions of our algorithm. 



4.4 Some illustrative results 

We illustrate our method on benzene, indigo and fullerene Ceo- Benzene serves as a test whether our basis represents 
the non interacting response function correctly. For indigo and fullerene, we compare our results with those of ADF 
[5] and Quantum Espresso |25j . respectively. 

In figure [4] we compare two theoretical calculations with experimental data for benzene (CeHg). Both calcu- 



lations make use of the iterative method of section 4.2 but we used (i) the basis of dominant products and (ii) the 
basis of products of (Kohn-Sham) eigenstates in these calculations. In the basis of products of eigenstates, the non 
interacting response xo is a diagonal matrix [27] of size N occ N v i T t ~ N 2 , that is easy to compute. Therefore, the 
comparison (i) vs (ii) allows us to assess whether there are enough dominant products for representing the Kohn-Sham 
response function. The size of the set of products of molecular orbitals grows as iV 2 and the TDDFT kernel becomes 
a dense matrix, leading to an overall complexity scaling 0(N 3 N LJ ) and limiting the possibility of such a comparison to 
relatively small molecules like benzene or naphthalene. The basis is chosen to contain 1238 dominant products, which 
is considerable less than the 108*109/2=5886 original products. 

The good agreement of the two theoretical calculations indicates the validity of our method as a whole and the 
adequacy of the dominant product's basis for representing the Kohn-Sham response function. 

By comparison with experiment, the collective peak is slightly blue-shifted. There are several reasons for this 
discrepancy. A first reason is the limited applicability of the simple LDA functionals — we use the Perdew-Zunger 
LDA functional, which is the default functional in the SIESTA package |21] . A preliminary calculation with a GGA 

1 In this equation and below, we focus on the calculation of a particular component of the polarizability tensor and remove tensor indices from 
the polarizability P(co) and Cartesian vector indices from the dipole momenta = f F^(r)ridr. However, a corresponding block-Lanczos 
algorithm exists and allows for a faster simultaneous calculation of all components of the polarizability tensor 
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Figure 4: The interacting polarizability of benzene computed with two iterative methods compared with the experimental 
spectrum of benzene. The calculations differ by their basis - dominant products versus products of eigenstates. The input 
of both calculations is generated in a SIESTA run. Although different types of product basis are employed, the agreement 
between the theoretical results is good. The experimental spectrum shows a slightly red-shifted collective peak. See the 
text for a discussion of possible reasons of this discrepancy. 



functional shows better results. A second reason is the incompleteness of the LCAO basis set — we employed a double- 
zeta polarized (DZP) basis which is the default basis in the SIESTA package. A third reason is our choice of the pseudo 
potential. We have been using the pseudo potentials published on the SIESTA Internet site |28| (LDA potentials of 
Troullier-Martins type adapted from ABINIT package [29])- Surely one can perform a fine tuning of the parameters, 
but this is beyond the scope of this paper. Therefore, the remaining calculations presented below are done with the 
same set of parameters. 
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Figure 5: Comparison of the absorption of indigo dye computed with the ADF package versus our method. DFT 
calculations are done independently in ADF and SIESTA. Nevertheless, theoretical spectra agree except for a factor of 
2 in the height of the HOMO-LUMO peak. In the range between 3-4 eV experiment and theory differ. The deviation 
might be caused by the presence of solvent in experimental setup. 

In figure [5] we compare two theoretical calculations with experimental data [32] for the indigo dye (C16N2O2H10). 
The first calculation is done within the ADF package [5 with parameters similar to SIESTA's default parameters, 
while the second calculation is done with our method (details of the calculation are collected above in the discussion 
of results for benzene). 

Both calculations agree except for the strength of the HOMO-LUMO transition. The experimental spectrum is in 
overall agreement with both calculations: three resonances are seen, while the middle resonance is probably disturbed 
by the presence of the solvent. 
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Figure 6: Theory versus experiment for absorption of fullerene Ceo- Experimental data are from |30) and we compare our 
result with those of [31] . The theoretical results agree with each other and disagree with experiment. This might indicate 
an inadequacy of the LDA exchange-correlation functional for large molecules or be due to the presence of a solvent in 
the experimental setup. Excitonic effects may be another reason for this discrepancy. 

We kept about 2100 dominant products and 256 frequency points in this calculation. The current implementation 
took 3.5 hours, if run on one thread, on an Intel Xeon processor at 2.50 GHz. 

In figure [6] we compare theoretical and experimental spectra for fullerene Ceo- The calculation by Rocca et al |31| 
agrees remarkably well with our calculation, while the experimental results [30] are blue-shifted compared with theory 
(details of the calculation are collected above in the discussion of results for benzene). 

The disagreement might be due to use of the LDA functional. A solvent usually gives rise to a uniform red shift of 
experimental data [33] which is not the case in this example. We believe that excitonic effects cause the discrepancy 
in this large molecule. 

We kept about 8700 dominant products and used 128 frequency points. The current implementation requires 18.4 
hours, if run on one thread, on an Intel Xeon processor at 2.50 GHz. 

The agreement between our calculations with those done by other authors and with experiment validates our 
construction of a basis in the space of products and our construction of the Kohn-Sham response function x<>- 

5 Hedin's GW self-energy in the product basis 

It is known that particle-hole interactions in organic semiconductors cannot be described adequately by TDDFT. Such 
systems, however, can be modelled by Hedin's GW approximation [3"4"1 13"5] . 

Let us show how to express Hedin's GW self-energy [36] in terms of our concepts {F x (r), V£ b }. Although the 
self-energy operator E should, in principle, be determined self-consistently using Hedin's set of integral equations, even 
a "one shot" approximation to the self-energy has been shown to improve the quasi particle energies compared to 
TDDFT energies. The self energy is given by [361137] 




where W = 1 _f xo /h is a RPA-screened Coulomb interaction that makes use of the non interacting KS response \o 
and Go is the non interacting Green's function. To see the form this equation takes in our basis of dominant functions, 
we define a self energy matrix and expand the Greens function Go 




Inserting £ = iGoW and using the representation of products (|3J) we find the following tensor form of the GW 
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approximation 

T, ab (t~t')=i V^ a 'v bb 'G a/b/ (t~t')W^(t^t'), 

[i,is,a' ,b' 




Figure 7: Diagram for the self-energy E in a "one shot" GW approximation. The vertices V^ a and V„ b are sparse - only 
the orbitals and products belonging to a given quadruplet of atoms contribute to the trace. Therefore, the computational 
effort for self energy E scales as 0(N 2 N LU ) once the screened Coulomb interaction is known. 

The direct computation of the screened Coulomb interaction requires a priori 0(N 3 N LJ ) operations while the 
rest of the calculation of E ai> scales as 0(N 2 N^). This is due to the locality of the vertex V^ b (see the diagrammatic 
representation of the self-energy in figure |7b. The situation here is similar as in the calculation of xJL- 



6 Conclusions 

In this paper we reviewed our extension of the LCAO method to densities and excited states and we gave first 
applications of this method. One application is a convenient construction of the non interacting Kohn-Sham response 
in 0(N 2 Nu) operations. Another application is the use of this response function to compute electronic excitation 
spectra within TDDFT linear response, again in 0(N 2 N tJ ) operations. We illustrated our method of computing 
spectra using benzene, indigo and fullerene and we also confirmed the 0{N 2 N LJ ) complexity scaling of our method. A 
drawback of our construction of xo is its high memory requirement — one could use a MPI parallelization to address 
this problem. We also mentioned that our method is suitable for application to the GW approximation which we cast 
in an appropriate tensor form. 

Because our method provides (i) a simple basis for the electronic density and (ii) the Kohn-Sham response function 
we believe it will be useful in treating excitonic effects in molecular physics. 
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